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Abstract 

This paper is concerned with the numerical approximation of the 
isothermal Euler equations for charged particles subject to the Lorentz 
force (the 'Euler-Lorentz' system). When the magnetic field is large, or 
equivalently, when the parameter s representing the non-dimensional 
ion cyclotron frequency tends to zero, the so-called drift-fluid (or gyro- 
fluid) approximation is obtained. In this limit, the parallel motion rel- 
ative to the magnetic field direction splits from perpendicular motion 
and is given implicitly by the constraint of zero total force along the 
magnetic field lines. In this paper, we provide a well-posed elliptic 
equation for the parallel velocity which in turn allows us to construct 
an Asymptotic-Preserving (AP) scheme for the Euler-Lorentz system. 
This scheme gives rise to both a consistent approximation of the Euler- 
Lorentz model when e is finite and a consistent approximation of the 
drift limit when e — > 0. Above all, it does not require any constraint 
on the space and time steps related to the small value of e. Numerical 
results are presented, which confirm the AP character of the scheme 
and its Asymptotic Stability. 
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1 Introduction 

This paper is concerned with the construction of a numerical scheme for the 
system of isothermal Euler equations for charged particles subject to the 
Lorentz force (which we will refer to as the Euler- Lorentz system). More 
precisely, we are interested in the regime where the inertia terms which bal- 
ance the pressure and Lorentz forces in the momentum balance equation are 
scaled by a small parameter e. The parameter e represents the inverse of the 
ion gyrofrequency around the magnetic field axis scaled by a characteristic 
time of the experiment. When e tends to zero, the so-called drift-fluid or 
gyro- fluid regime is reached [211 1211 • 

In the drift-fluid approximation, particles are confined along the magnetic 
field lines. As a consequence, the dynamics along the magnetic field lines is 
much quicker than across it. In the limit e — > 0, the parallel motion assumes 
an instantaneous equilibrium in which the pressure force equilibrates the elec- 
tric force. This equilibrium is attained through acoustic waves propagating 
at infinite velocity in a similar fashion as what happens in a low Mach num- 
ber fluid. These acoustic waves adjust the parallel velocity instantaneously 
in such a way that this equilibrium is satisfied at all times. In this way, the 
parallel velocity plays the role of a Lagrange multiplier of the constraint of 
zero total aligned force. The first goal of this paper is to give an equivalent 
formulation of the drift-fluid approximation that enables us to calculate this 
parallel velocity (contrary to what is sometimes written in the literature [2TJ , 
it is possible to find such an equation). 

The second goal is to design a valid scheme for both regimes e ~ 1 and 
e — > 0. This scheme gives rise to both a consistent approximation of the 
Euler-Lorentz model when e is finite and a consistent approximation of the 
drift-fluid limit when e — > 0. Above all, it does not require any constraint 
on the space and time steps related to the small value of e. This type of 
schemes is usually referred to as Asymptotic Preserving schemes (AP). 
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Asymptotic Preserving schemes have been proposed in a variety of con- 
texts, such as hydrodynamic or diffusion limits of kinetic model [21 EH1 EIH E21 
|3T| HI [T7], relaxation limits of hyperbolic models [23j EH [18], relaxation lim- 
its of Complex-Ginzburg-Landau equations (TO], low-Mach number limits of 
compressible fluid models [9]. In the plasma physics context, these schemes 
have appeared in relation with the quasi-neutral limit of the Euler-Poisson 
system [HI ELI] or Vlasov-Poisson system [HIE]- 

Such schemes are of great potential interest to the simulation of strongly 
magnetized plasmas such as those encountered in space plasmas or in Toka- 
mak devices like ITER. First of all, there are several advantages to using the 
original Euler-Lorentz model instead of the limit drift-fluid model. Indeed, 
the drift-fluid model is a mathematically complex model: the constraint of 
zero total force makes it a mixed-type model, with certain characteristics of 
an elliptic problem, like the incompressible Navier-Stokes equation. Dealing 
with this constraint is numerically challenging, and is at least as difficult as 
dealing with the incompressibility constraint in the Navier-Stokes equation. 
In the literature, various drift-fluid models have been proposed on physical 
grounds [H [21 [HI HSl HEl UHl I2DI [2HI EE EOl E21 E31 EH]- However, their rela- 
tions to the drift-fluid model which can be derived from the formal asymptotic 
analysis developped below are unclear. This is why we find preferable to rely 
on the original Euler-Lorentz model, in which the momentum conservation 
equation directly follows from first physical principles. 

Another advantage of AP schemes is their ability to deal equally well 
with the asymptotic regime e — > and the 'normal' situation e = 0(1). This 
is potentially very interesting for situations in which part of the simulation 
domain reaches the asymptotic regime and part of it does not. The usual 
approach for dealing with such occurences is through multiphysics domain 
decomposition: the full Euler-Lorentz model is used in the region where 
e = 0(1) and the drift-fluid limit model is used where £ < 1 (we assume 
the dimensionless parameter e is computed using local estimates of the mag- 
netic field strength and can be considered as a function of the space and 
time coordinates). There are several drawbacks in using this approach. The 
first one is the choice of the position of the interface (or cross-talk region), 
which can influence the outcome of the simulation. If the interface evolves in 
time, an algorithm for interface motion has to be devised and some adaptive 
remeshing has to be implemented, which requires heavy code developments 
and can be quite CPU time consuming. Determining the right coupling strat- 
egy between the two models can also be quite challenging and the outcome 
of the numerical simulations may also depend on that choice. Because these 
questions do not have a straigthforward answer, multiphysics domain de- 
composition strategies often lack robustness and reliability. Here, using the 
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original model with an AP discretization method everywhere prevents from 
having to introduce questionable physical artefacts in the model and permits 
to use the same code everywhere in both regimes. 

The potential application of these methods is the simulation of fluid tur- 
bulence in Tokamak plasmas. There has been considerable literature on this 
problem [TJ [21 EH HH EH [213 ES] • The present work is far from being at a com- 
parable development stage, since the present numerical tests are restricted 
to given uniform magnetic and electric fields in a two-dimensional setting. 
Nonetheless, this is an unavoidable intermediate step to check the perfor- 
mances of the method. The numerical tests being successful!, this approach 
will soon be extended to an arbitrary given magnetic field and coupled to 
the dynamics of the electron fluid through quasineutrality assumptions. 

The assumption that the ion fluid is isothermal is only made for simplicity. 
An energy equation for the ion fluid can be considered instead. The approach 
extends easily to this case. We will report on it in future work. 

This paper is then organized as follows. The isothermal Euler-Lorentz 
model in the drift-fluid scaling and the drift-fluid limit are presented in sec- 
tion [21 An Asymptotic Preserving time discretization of the isothermal Euler- 
Lorentz model in the drift-fluid scaling and a full discretization of this scheme 
in a reduced 2-dimensional setting are proposed in section [3] Numerical re- 
sults are given in section H] and finally, conclusions are given in section [51 

2 The isothermal Euler-Lorentz model and 
the drift-fluid limit 



2.1 The Euler-Lorentz model 

We are concerned with the Euler-Lorentz model describing the isothermal 
flow of positive ions in a tokamak. In this model, we neglect the electrons 
and suppose that the electric and magnetic fields are given. In the drift- 
fluid asymptotics, we let e be a typical scaled value of the gyro-period of the 
particles, i.e. the period of their rotation motion about the magnetic field 
axis. The scaled isothermal Euler-Lorentz model takes the form: 



d t n £ + V • (n £ u £ ) 







d t (n £ u £ ) 



V • (n £ u £ ® u £ 



(2.1) 

TVn £ =n £ (E + u £ x B) , (2.2) 



where n £ , u £ are the density and the velocity of ions, respectively. The 
quantity T is the ion temperature. Here, the electric field E and the magnetic 
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field B are assumed to be given functions. The symbol V is the gradient 
operator while V- denotes the divergent operator. 

This scaled model is obtained from the unsealed Euler-Lorentz model by 
introducing characteristic scales for length xq, time to, velocity uq, density 
n , temperature T , electric field E and magnetic field B . As usual, we set 
Xq = Uot and the characteristic electric and magnetic fields are assumed to 
follow the relation E = u B , so that the gyro-frequency of the ions is given 
by uj = qBo/m = qEo/(muo) (where q is the ion electric charge). In doing so, 
two dimensionless parameters appear, the Mach number M = uq/c s where 
c s = (To/m) 1 / 2 is the sound speed (and m is the ion mass) on the one hand, 
and the scaled gyro-period e = mj (qB t ). In the drift-fluid asymptotics, we 
assume that the Mach number and the gyro-period are linked by M = y/e, 
which leads to the scaled problem (12.1 1) . (12.21) . 

The following notations will be useful: the director of the magnetic field 
is denoted by b = B/B where B is the Euclidean norm of B. 

Any vector quantity v can be split into its parallel (||) and perpendicu- 
lar (_L) parts as follows: 

v = v\\ + v± = v nb + v± , v\\ = v ■ b , v±_ = b x [v x b) . 

Next, we introduce the parallel gradient Vy0 = b ■ V0 for any scalar func- 
tion 0. The quantity Vy0 is a scalar. In the same way, we also introduce the 
parallel divergence, given for any vector field v by Vy ■ (v\\) = V • (v\\b). This 
operator can be related to the parallel gradient by following equalities, 

V||^j|=V.(|i?)=BV||(|), (2.3) 

since the magnetic field is a divergence-free vector. We can write V ■ v = V ■ 
u_i_+V|| ■ (f||) . Note that Vy • (vy) is also a scalar. More generally, we introduce 
the parallel divergence Vy ■ <fr, acting on a scalar <p, by Vy • <f) = V • (<frb). The 
operators Vy and Vy are formal adjoints. Let us consider two scalar- valued 
fonctions <fi and ip defined on a regular domain Q, and let us assume also that 
(f) and if) vanish on the boundary dQ, for simplicity. We have, 

J V ll( V l|-^) <l>dx = J b ■ V ( V ■ W ) i> dx 

V • (06) ) (V • (#) ) dx 

J Jvwt) (vir^) dx. 
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2.2 The Drift-fluid Limit 

The formal limit e — > in the isothermal Euler-Lorentz model ( 12. II) . ( 12.21) . 
leads to the so-called isothermal drift-fluid model: 

d t n + V • (nti) = , (2.4) 
T Vn = n(E + u x B) . (2.5) 

The constraint (12.51) completely determines the velocity u. Indeed, taking 
the parallel and perpendicular components of (12.51) leads to 

nu ± = \-b x (T Vn - nE) , (2.6) 
B 

TV\\n-nE\\ = 0. (2.7) 

After dividing by n, we find that the first term at the right-hand side of ( 12.61) 
is the diamagnetic drift velocity while the second one is the E x B drift 
velocity. 

Eq. (12. 7p can be recast in the form of an elliptic equation for u\\. Indeed, 
(12.41) can be written: 

) + V|| • (nu\\) = 0, (2.8) 

Applying Vy to (12.81) . noting that [Vy, d t ] = —d t b- V (where [-, •] denotes the 
commutator) and inserting ( 12. 7p leads to: 

- Vy ( Vy • (n«||)) = dt - d t b ■ Vn + Vy ( V ± • K)) . (2.9) 

This is a one-dimensional elliptic equation for u\\ along the magnetic field 
lines, which is well-posed through the above mentioned duality between the 
parallel gradient and parallel divergence operators. Therefore the parallel 
component u\\ can be computed explicitly through the resolution of this el- 
liptic equation, provided that adequate boundary conditions are given. The 
boundary conditions depend on the specific test case under consideration. 
They will be discussed in the numerical section below. 

The drift-fluid model consists of equations (12.41) . (12.61) and 



2.3 A reformulation of the isothermal Euler-Lorentz 
model 

The scaled Euler-Lorentz model in the drift-fluid asymptotics is a singularly 
perturbed problem: in the drift-fluid limit, the type of certain equations 
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changes. Indeed, in the Euler-Lorentz model the velocity is given by a time 
evolution equation of hyperbolic type while in the drift-fluid limit, the per- 
pendicular velocity is given by an algebraic equation, while the parallel com- 
ponent is found through solving an elliptic type equation. To find an AP 
scheme, it is essential to 'regularize' the perturbation, i.e. to reformulate the 
Euler-Lorentz model in such a way that the limit equations for the velocity 
appear explicitly in the system of equations. The goal of this section is to 
find such a reformulation. 

For the perpendicular component of the momentum, we take the cross- 
product of ( 12.2ft with b, which leads to: 



B (n e u e ) 



edAbx (n F u F 



b x 



+ e 



dtb) x (n F u F 



T Vn £ + n £ E 
b x ( V • (n F u F 



(2.10) 

Formally, when e — > in eq. (I2.10p . we recover the equation for the perpen- 
dicular component of the momentum in the drift-fluid limit model (12. 6|) . 
We now take the scalar product of (12.21) with b: 



0, 



d t b) ■ (n £ u £ ) + b- ( V ■ (n £ u £ ® u 



T Vn £ + n £ E 



'XIV, 



Since u\\ cannot be computed explicitly from this equation in the limit e — > 
0, we are led to reformulate the equation ( 12.111) . We first take the time 
derivative of ( 12. lip and get 
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\n £ u £ 



u F 



I - d t [[dtb) ■ (n £ u £ )J +d t [b- [V 

= d t {n £ E\\) - Td t b ■ Vn £ -TV\\d t n £ . (2.12) 

Now, applying Vn to (12.11) (rewritten in the same fashion as (12. 8p ) leads to 



edt 



[n £ u £ ) 



TV||(V|| 



[n F u 



d t b' 
dt{n £ E, 



,n £ u £ 



ed t [b- V • (n £ u £ ® u £ 



Td t b ■ Vn £ + TV||(V_l ■ (n £ u £ 



(2.13) 



We notice that eq. (12. 9p is the formal limit of eq. (I2.13P when e — > 0. 
Eq. (I2.13P is a wave equation for u\\ associated with the elliptic operator 
Vy ■ (V||), which is well-posed provided that suitable boundary conditions are 
given. This wave equation describes the propagation of disturbances along 
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the magnetic field lines, which propagate at a velocity of order 0(e^ 1 ^ 2 ). 
In the limit e — > 0, an equilibrium described by (12.91) is instantaneously 
reached through waves propagating at infinite speed. Eq. ( I2.13P provides an 
equivalent formulation to ( 12.2ft for u\\, but which does not become singular 
when e — ► 0. 

Therefore, the reformulation of the Euler-Lorentz model consists of eqs. 

([21]), (MD, 0233D. 

3 An Asymptotic Preserving scheme for the 
isothermal Euler-Lorentz model in the drift- 
fluid approximation 

3.1 Time semi-discrete scheme 

The purpose of this section is to build an AP scheme for the Euler-Lorentz 
model, i.e. a scheme which is consistent with the Euler-Lorentz model when 
e = 0(1) and with the drift-fluid limit model when e <C 1. The AP property 
mostly relies on an appropriate time-discretization. We will investigate this 
point first. Of course, we have in mind that time semi-discrete schemes of 
hyperbolic problems are unstable unless some diffusion is added. In this 
section, we assume that the gradient operators are actually approximate 
operators which encompass the requested numerical diffusion. The space 
discretization is discussed in detail in section 13.21 

Our AP time semi-discrete scheme relies on use of the reformulated equa- 
tions ( 12.11) . ( 12.1 Oft . ( I2.13p . However, rather than looking for a discretization of 
them, it is more efficient to start from a discretization of the original formu- 
lation (12.11) . (12. 2p and find a scheme which allows the same reformulation as 
the continuous problem and the derivation of the discrete equivalent to (12.11) . 
(12.1 Op . (12.131) . In this way, we are guaranteed to find a suitable discretization 
also in the regime e = 0(1) which we could miss otherwise. 

We first introduce some notations. Let B m be the magnetic field at 
time t m , B m its magnitude and b m = B m / B m its director. For a given 
vector field v, we denote by (v )™ and (v)™ its parallel and perpendicular 
components with respect to b m . Similarly, we denote by V™ and V™- the 
parallel gradient and divergence operators respective to this field. 

To calculate the solution of the isothermal Euler-Lorentz model in the 
drift-fluid approximation ( 12.11) . ( 12. 2ft . we propose the following time semi- 
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discrete scheme, 



„m+l _ „m 

e + V • (n £ u £ ) m+1 = 



\ m 



At 

{n £ u £ ) m+1 - (n £ u £j 

h V • (n £ u £ g> n £ ) 



At 



+ T ( Vnf ) 



m+l 



ii"£ m+1 + (n £ u £ ) m+1 x B m+1 



Here, the quantity (Vn*) m+1 is given by, 

(Vn#) m+1 = (Vn™)™ +1 + (V< +1 )7 +1 



(3.1; 



(3.2) 



(3.3) 



In this scheme, the mass flux, the parallel component of the pressure force and 
the Lorentz force are evaluated implicitly while the perpendicular component 
of the pressure force is evaluated explicitly. We show that these choices 
permit a reformulation of the scheme into discrete equivalents to eqs. (12.11) . 

(EIUD, (J2D3D- 

We first investigate the transverse component and take the cross-product 
of (HOD with b m+1 . This leads to, 



[n £ u £ 



-l £ 



'771+1 



x m + 1 



At B m+1 



'm+l 



— \ n £ u £ ) 
I At 



eV ■ (n £ u £ ® u £ ) rn - T Vn™ + n™E m+1 



(3.4) 

which is a discretization of eq. (12.101) . where (dtb) x (n e ii e ) ~ ((b m+1 — 
b m )/At) x (n £ ™ £ ) m 

We now compute the scalar product of (13.21) with b m+l . We get: 



^((n £ ^r +1 )7 +1 +TV 

e 



m+l m+l 

II £ 



■m+l 



At 



(n £ u £ ) m - e[ V ■ (n e n £ ® ™ £ ) m ) + <^ m+1 



(3.5) 



which is a discrete version of eq. (12.111) . Differentiation of the discrete mass 
conservation equation (13.11) in the parallel direction gives 



^m+l^m+l _ ^m+l^m 



(3.6) 



Atv- +i (v-((^ £ r +i )™ +1 ) 

- AtV|p +1 (vjp +1 ■ ((n e u e ) m+1 )J +1 ) , 
which can be used to eliminate n™ +1 in favor of {n £ u £ )l +l in (ESD. This 
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leads to: 



At 



((n £ u £ ) m+1 )™ +l -T AtV™ +1 (v™ +1 ■ ({n e u e ) m+1 )™ +1 ) 



TAtVf +1 (y ■ {(n £ u £ ) m+1 ) 



m+l 



— (n £ u £ ) 
At y 1 



el V • {n £ u £ ® u £ ) 



TV™ +1 n m 



(3.7) 

m+l 



This equation is a one dimensional elliptic equation (along the magnetic 
field lines) for the quantity ((n £ u £ ) m+l )™ +1 . It is the discrete counterpart 
of fl2~T3D but the link with (123311 is not fully direct. Eq. (J33D is rather a 
discretization of the following equation: 



ed t (n £ u £ ) 











TV,, 


[/V 


((™e«e)|| 













e[(^6) • (n £ w e ) - 6- (V (n £ u £ ®u 



TVn 



V ■ (n e u £ ) 



b- 



T Vn™ - n F E 



ds 



(3.1 



which is obtained through the reformulation process outlined in section 12.31 
when the mass conservation equation is used in time-integrated form 



n £ = n™ 



f V ■ ((n e u e )j_) ds - f V,| • ds. (3.9) 

Jt m Jt m 



That (13.81) is equivalent to (12.131) is easy and is left to the reader. 

Now, we investigate the limit e — > in ( 13.41) . ( 13. 7ft . leaving At unchanged. 
We get, 



[nu 



\m+l 



m+l 



X 



(3.10) 

+ iym+l 



T Vn m + n m E m+1 

T At V™ +1 (y™ +1 ■ ((nu) m+1 )™ +1 ) = T At ■ ((nu) m+1 ) 

-TV^ +1 n rn + [n m E m+1 }™ +1 . (3.11 



This is the discrete counterpart of the drift-fluid equations ( 12.61) . ( 12.91) . 
Therefore, the limit e — > can be taken in the scheme ( 13. ip . ( 13.41) . ( 13. 7p and 
the resulting scheme is consistent with the drift-fluid equations. This shows 
that the time semi-discrete scheme (13.11) . (13.41) . (13.71) provides an Asymp- 
totic Preserving discretization of the Euler-Lorentz model in the drift-fluid 
limit. This scheme enables us to compute the solution of the isothermal 
Euler-Lorentz model for all regimes ranging from e = 0(1) to e <C 1 with 
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the same time step At. A conventional scheme would require to let At — > 
simultaneously with e — > 0. An AP scheme is free from this constraint. 

As a comparison, let us investigate the conventional time semi-discrete 
scheme for the Euler-Lorentz model which is the following: 

^m+l „m 

~V + V .(n £ u £ r = 0, (3.12) 



At 

m+1 



[n £ u £ - [n £ u £ ffl ,™ 
h V ■ [n £ u £ ® u £ ) 



At 

^m+l^m+l + ( rleW£ )-+ 1 x . (3.13) 



e 



The difference with our scheme is that the mass flux and the pressure force 
are both evaluated explicitly. We note that the Lorentz force at the right- 
hand side of (I3.13P is still implicit otherwise some obvious instabilities arise 
in the discretization of the cyclotron rotation (the ux B term in the Lorentz 
force). 

In taking the formal limit e — > in this scheme, we find that the parallel 
component of the velocity is now defined by the constraint 

TVf +1 n m - (n m+1 E m+1 )™ +1 = . (3.14) 

We cannot reproduce the computations leading to the elliptic equation for the 
parallel velocity (13.111) (or an analogous equation), because (13.14j) provides 
the value of V™ +1 n m instead of that of Vjp +1 n m+1 . So, we are bound to 
using the discrete time derivative at the left hand-side of (13.131) to update 
the value of {nu)\\ which requires At < e. Therefore, the conventional time 
semi-discrete scheme is not Asymptotic-Preserving. 

We also see that in our AP scheme, we need to evaluate both the mass flux 
and the parallel component of the pressure force to get an elliptic equation 
for (nu) || . If the pressure force alone is taken implicitely, this results in an 
equation for (nu)n which is ill-posed. 



3.2 Fully discrete scheme 

A two-dimensional case. For the sake of simplicity, we restrict ourselves to 
a 2-dimensional case with a constant in time, uniform in space magnetic field 
lying in the computational plane. Therefore, we assume that the magnetic 
field is directed along the y-axis and that the plasma lies in x, y-plane, with 
translation invariance in the ^-direction. However, a possible non-zero plasma 
velocity is assumed in the z direction. In these conditions, the isothermal 
Euler-Lorentz model in the drift-fluid asymptotics (12.11) . (12. 2p is written (we 



11 



now omit the indices e for the sake of simplicity) , 

d t n + d x {nu x ) + dy{nuy) =0, 

+ T d x n = nE x — nu z B y , 



\nu x u y ) 



T d y n = nEy , 



d t {nu x ) + d x (nu 2 x ) + d y ( 
d t {nuy) + d x {nu x u y ) + d y (nu 2 y ) 
d t (nu z ) + d x (nu x u z ) + d y (nu y u z 

Then the time semi-discrete AP scheme for the system (I3.15P reads 



(3.15) 



nE z + nu x B y , 



[nu n 



\ m+l 



£ 1 

AtB 



6 1 

At 5 v 



(n«x) m+1 + (mt 2 ) 



\ m+l 



m+l 



1 

1 



At 



e[d x (nu x u z ) 



+dy(nu y u z ) 



n m E ? 



At 



rat. 



2\ m 



+d y (nu x u y 



d x (nu 2 x ) 



Td x n m -n m E x 



(3.16) 



(3-17) 



- { nu y 



i m+l 



TAtdy (dy (nUy) m+1 ) = TAtdy {d^UU^) 



LAt [nUy} 



2\m 



d x (nu x u y ) m + dy(nu 2 y ) 



Td y n r 



+ d x (nu x ) m+1 + dy(nu y ) 



At 



m+l 







+ n m ^3.18) 
(3.19) 



where B = b ■ B = B y . We now give the full space-time discretization based 
on the above discussed AP scheme for system f)3.15p . 

Fully discrete scheme in the two-dimensional case. For numeri- 
cal purpose, let us consider a Cartesian mesh of the calculation domain 
(zi-i/2, x i+1 / 2 ) x (yj-i/2, yj+1/2), i,j = 1-N. Then eqs. (|3.16|), (|3.17|), 
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(I3.18P and (13.191) are discretized according to, 



(nu x ) 



m+l 



£ 1 



raz., 

At B 



m+l 



+ (rat.. 




1 

(nUrU 



—— (nu z ).- 
At v Aj 



nu y u z)ij +1 / 2 



\ m+l 



1 



2 T 

to x H n 



i+l/2j 



zJi-1/2 j 



{nUyU z J i j_i/ 2 



■-— \nu x ).. 
At v JlJ 



2 T 

nu x H n 



(n m E z ) tJ 



(nu x u y )™ +l/2 - (n« x « w )™._ 1/2 ) - (n m £ x ) 



£ / \m+l 

At v yAj 



TAt (aJ(n« w ) ra+1 ) i . = TAt (a,(a 



£ / \m ^ ( ( \Tn / \m \ 

At ^ nUy '^ ~ ~Ax \ ( nUxU yn+i/2 j ~ y nu - u y)i-i/2j ) 




m 



mt. 



m + l 



)) 



'J 



r 



n«„ H « 



/ \m+l 
\ nu x)i-\/2j 



(3.20) 



(3.21; 



(^)J +1 (3.22) 



m+l 
VJij+l/2 



\m+l 
VJij-1/2 



(3.23) 



Here n^ +1 is the density in the cell (xj_i/ 2 , x i+ i/ 2 ) x (yj-i/2, yj+1/2) at 
the time t m+1 . The quantity ((nu x )™ +1 , (^w.z)^ +1 ) is the perpendicular 
part of the momentum while is the parallel part of the momentum 



in the cell (x. 



i-l/2j ^i+1/2 



.) x (yj_i/2, yj+1/2) at the time t m+1 . The terms 
(")i+i/2j arL( ^ (0^+1/2 denote the numerical fluxes at the interfaces Xi+1/2 and 
Vj+i/2 of the corresponding quantities, respectively. The second order terms 



(dy(nu y ) m+ ) and (d y (d : 



\ m+l 



) ) will be discussed below. 



Discretization of the hyperbolic part. To calculate the numerical fluxes 



at the interfaces 



h+i/2j 



and 



we use the Pq scheme [12]. To be 



more precise, let us consider the interface £1+1/2 separating data U™ 1 , Z^+ij 
for the corresponding Riemann problem at time t m , where 
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Let us denote by M^^j = (j^i+i/2j^ nu T+i/2j^ the average state between 

the states U^, Ui+i j- The average state U^^ - is the Roe average state [M] 
given here by following formula: 



7"< 



ivi ' ' " — /in in 777 

n i+i/2j - v 'i Wi' 



i+1/2j ,frE+,/W 



and the momentum of the average state is reconstructed as 

771 ^~77l ^771 

nU i+l/2j — n i+l/2j U i+l/2j ■ 

Then the numerical fluxes are given by, 

[nu x u z ) ij + [nu x u z ) i+lj 



{nu x u z ) i+ y 2 j 



f (nu z )T +ll -(nu z )-) ,(3.24) 



nu x H n 



2 

T \ m ( T 

nut. H n + nut H n 



rp \ rri \ F I V F 

2 , 1 \ \ t /ij\ c / i+lj 



i+l/2j 



m 

a i+l/2j 



{nu x u y ) i:j + (nua-tty) 



2 



- ,(3.25) 



nuy)^ - (™%)™J ,(3-26) 



\m+l . / \m+l rn 



(nu x )T^ /2l = " ' l3 2 " - (<U; - n$ , (3.27) 

/ \m+l i / \m+l m 

nuj^ , = - * J - " - fnS., , - n£) . (3.28) 



"yJi+i/2j 2 2 v 

Here, the speed a™!/ 2 j is given by 



m / i m,— i i m,+ i \ / on \ 

Oi+i/2,- = max |a ' . 1 |a. ' 1/2 J . (3.29) 



b i+l/2j — L1LaA \ \ u i+l/2j\i ^1+1/2 j 

with 

771. — • / rn p ^rrt, F 

°i+l/2,- = mln ~ C ' M * i+1/2J ~ C 

a^ /2j =maxfe +i/2j + C ^,^r +1 ,+c £ 



(3.30) 
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and c £ = ■\frjs standing for the sound speed. The numerical fluxes across 
interfaces j/j+1/2 are computed similiarly. 



Discretization of the second order terms. The second order terms are 
computed with centered spatial discretizations: 



niiJ..,, —zlriUv)-. +[nu v ).. , 
\y y \nuy) ).. - ^ 2 , 



and 



(d y (d x (nu x ) m+ )).. = — [ (d x {nu x ) m+ ). . +1/2 - (d x (nu x ) m+ ).._ y 
where 

d x ( nu y) m+1 ). j+1/2 = 2 ( d * ( nu *) m+1 ) ij+1 + 2 ( n ^) m+1 )i^ 

m+l / \m+l 



1 (™«*)™i J+ l - (™.) 



t-lj'+l 



+ 



2 2 Ax 

/ \ra+l / \m+l 

1 {nu x ) i+lj - [nu x ) i _ lj 



2 2 Ax 



To solve the elliptic equation for {nu y ), suitable boundary conditions need 
to be specified. In particular, using above discretizations eq. (I3.22p can be 
recast in the form 



A x rn+1 = nns m+1 , (3.31) 

where A = A (e, At, Ay, T) is a regular matrix, X m+l = ((m^)™ +1 ) ij=1 jy 
is the vector of the parallel momenta in all the computational domain and 
nnS m+i = nnS m+1 (e,At,Ax,Ay,T,(nu x ) m ,(nu y ) m ,E™) is the right- 
hand side, which is known. The linear system (13.311) is solved by a Gaussian 
elimination method with partial pivoting. 

Choice of the time-step. As usual, the time-step At is chosen such that 
the CFL condition is satisfied. In 2D geometry, this condition takes the 
following form, 

At ( max C 1/2 , + -^ max a"! +1/2 J = CFL < 1 . (3.32) 

\Ax i<i,j<N l+1/2] Ay i<i,j<N % i+ v l 1 ) v ' 

where a™ 1/2j is defined by (l3~2^D . The scheme f[3~2U]) - ([3~23l with the time- 
step algorithm (13.321) will be referred to as the resolved AP scheme. 
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When e tends to 0, the sound speed c £ = y/T/e takes large values and 
the time-step calculated from (13.321) becomes very small. Due to the AP 
character of the scheme, we do not need to constrain the time step to stay of 
the order of e. Therefore, the sound speed can be removed from the definition 
of the velocities used to compute the numerical viscosity in the interfacial 
fluxes. These new velocities are defined as follows: 

j = maX (l 5 m/2jl> l°i+l/2jl) ' (3.33) 

with 

~m,— ■ I m -~m \ ~m,+ 

a i+X/2j = mm ,. J > °i+l/2j 

and a new CFL condition is introduced: 

At (-!- max aS. 1/2 , + -!- max a"! +1/2 ] = CFL < 1 . (3.35) 

\Ax l<i,j<N l+1/2] Ay l<i,j<N J - y ' 

It is important to notice that, with this new expression, the time step can 
be chosen independent of e. 

The scheme f l3.20l) -( l3.23l) where velocities a™ l j 2 - are substituted by a™ x , 2 
given by eq. f)3.34p and with the time-step algorithm (13.351) will be called the 
non-resolved AP scheme. The numerical tests will show that this choice gives 
rise to a correct solution. 



= max(^ +i/2 .,<^.J , 

(3.34) 



4 Numerical tests 

4.1 Geometry and test case 

Our test-case is two-dimensional: the physical quantities depend only on the 
coordinates x, y. The magnetic field is assumed to be uniform and directed 
along the y-axis i.e. B = (0, B y , 0) while the electric field is directed along 
the 2;-axis, E = (0, 0, E z ). The test-case geometry is depicted in Fig. [TJ 

The domain is a square of side 1 as shown in Fig. [2] while the parameters 
are given in the table [TJ the initial values of physical quantities are put in 
the Q column while the boundary conditions are given in the columns /, 77, 
/// and IV, accordingly. 

For the considered test case, the exact drift-fluid approximation is sta- 
tionary and uniform in the whole domain, and given by 

n = 1 , nu x = —1 , nu y = 1 , nu z = , T = 1 (4.1) 
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Figure 2: Schematics of the computational domain. 



We observe that the initial and boundary data are 'well prepared': they are 
perturbations of order e of the drift-fluid limit. Indeed, if 'unprepared' initial 
and boundary data are used, large (of order 1) initial and boundary layers 
appear in which the exact solution is significantly different from the drift-fluid 
limit. In order to correctly capture these initial or boundary layers, there is no 
other way than using a time and space resolved scheme in which both the time 
and space steps are of order e. However, the goal of an AP scheme is not to 
capture the initial and boundary layers accurately, but to provide a consistent 
approximation of the correct drift-fluid limit where it applies, i.e. away from 
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Table 1: Simulation parameters for the test-case. 





Vt 
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J J 


III 


IV 
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1 


1 + e 
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1+e 
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nu x 





-1 
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-1+e 


-1+e 


nu y 
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1 + e 


1+e 


1 


nu z 








cr 
c 





e 


T 
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1 


1 


1 


1 


By 


1 


1 


1 


1 


1 


E z 


1 


1 


1 


1 


1 



these layers. An accurate verification of this property requires a test solution 
which is not polluted by the initial and boundary layers and therefore, the 
need of well-prepared initial and boundary conditions. In section 14.41 for 
the sake of completeness, we show some numerical results with unprepared 
initial and boundary conditions. 

4.2 Simulations for e <C 1 

Here we would like to demonstrate that the AP scheme is consistent with 
the drift-fluid limit even for large time steps compared to e. By contrast, 
the conventional scheme is shown to be unstable for time steps larger than 
e. Three values of the parameter e are used: e = 10 -5 , e = 10 -6 and 
e = \/e mac hine = 1.5 10~ 8 . For these values, we observe the numerical solu- 
tion at times 1, 0.1 and 0.01, respectively. The reason for choosing smaller 
observation times when e is smaller is due to the increase of computation 
time when the time step resolves e. The CFL number is taken equal to 0.5 
for all simulations. A uniform mesh is used for both the x and y direction 
with steps Ax = Ay = 0.01. 

Resolved case. We first compare the conventional and AP schemes in the 
resolved case, i.e. when At is smaller than e. The results given by both 
the conventional and AP schemes are displayed in Figs. [3] for e = 10~ 6 and 
compared with the drift-fluid limit. Here, as well as in the forthcoming 
pictures, the various physical quantities (i.e. the density n and the three 
components of the momentum nu) are shown as functions of x for a given 
value of y = 0.5. The computed solutions are indistinguishable and very 
close to the drift-fluid limit. 

In this case £ < 1, we want to test the consitency of the scheme with the 
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drift-fluid limit. To do so, we compute the difference between the numerical 
solution and the analytical solution 04. ip . This is not the numerical error in 
the conventional sense, since we do not compare the solution with the exact 
solution for the same value of e but rather with the exact solution in the limit 
e — ► 0. For this reason, we do not call this quantity, 'the error', but rather 
'the difference with the drift-fluid limit'. We normalize this difference by the 
exact value of the drift-fluid limit, except for nu z since this value is exactly 
zero. We display these quantities as functions of (x, y) for the resolved AP 
scheme with the value e = 10~ 6 on Fig. HI The picture is almost the same 
if we replace the resolved AP scheme by the resolved conventional scheme. 
For this reason, the latter is omitted. 

The maximal relative difference between the computed solution and the 
exact drift-fluid limit on the density and momenta for three values of e are 
given in table [2] for the resolved conventional scheme and in table [3] for the 
resolved AP scheme. They both show a very good agreement with the drift- 
fluid limit, which increases as e decreases. This is an expected result since 
both the resolved conventional and AP schemes use At = 0(e). For this 
range of time-steps the resolved conventional scheme is only slightly more 
accurate than the resolved AP scheme. 

Unresolved case. We now examine the unresolved situation, where e <C 
At. Results obtained by the non-resolved conventional and AP schemes for 
e = 10~ 6 are displayed in Fig. [5j We recall that "non-resolved" means that 
the viscosities are computed through (13.331) and the time-step through (13. 35ft 
instead of (13.291) . (I3.32p in the resolved case. 

Clearly, the computed solutions with the non-resolved conventional scheme 
are unstable, while those calculated with the non-resolved AP scheme remain 
stable and consistent with the exact drift-fluid limit. The difference between 
the computed solution and the exact drift-fluid limit on the density and mo- 
menta (relative difference for n, nu x and nu y and absolute difference for nu z ) 
of the solution for e = 10~ 6 are given in Fig. [HI for the non-resolved conven- 
tional scheme and in Fig. [7] for the non-resolved AP scheme. They confirm 
the stability and accuracy of the non-resolved AP scheme and the instability 
of the non-resolved conventional scheme. 

In Table HI the difference between the computed solution and the exact 
drift-fluid limit for the non-resolved conventional scheme for the three values 
of e are given, and similarly for the non-resolved AP scheme in Table El 
Again, the consistency of the non-resolved AP scheme with the drift-fluid 
limit on the one hand, and the instability of the non-resolved conventional 
scheme on the other hand, are confirmed. 

Fig. [HI shows the evolution of the time step with respect to time in the 
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case e = 10 -6 for both the unresolved and resolved AP schemes. We recall 
that the time step is not fixed once for all, but is recomputed at each time 
step using the CFL condition. However, we can see that, in log scale, the time 
step remains about constant. The time step for the unresolved AP scheme 
is about 3 decades larger than for the resolved AP scheme. In Table El we 
compare the time-step to the scaled gyro-period and we notice that it is of 
the same order of magnitude for the resolved AP scheme, as it should, and 
it is much larger (up to 4 decades !) for the unresolved AP scheme. 

It is even more interesting to compare the CPU time between the un- 
resolved AP scheme and the conventional scheme. Indeed, the AP scheme 
involves more complex computations than the conventional scheme, such as 
the inversion of linear systems. It is therefore legitimate to wonder whether 
this additional work does not completely balance the gain obtained through 
the use of larger time steps. To check this point, the CPU times for the 
conventional and non-resolved AP schemes for three values of e are given in 
table We see that the gain in CPU time is up to almost 4 decades with 
the smallest value of e. The gain scales about like \fe as it should. 

These comparisons show that the proposed AP schemes are very powerful 
in handling the numerical approximation of drift-fluid asymptotics. 

Table 2: Maximum of relative difference between the computed solution 
and the exact drift-fluid limit (%) from the resolved conventional scheme 
on the density n, x-component of the momentum nu x , y-component of the 
momentum nu y and absolute difference between the computed solution and 
the exact drift-fluid limit on the z-component of the momentum nu z . 



e 


n 


nu x 


nu y 


nu z 


1(T 5 


0.0087 


0.00714 


0.145 


0.0174 


1(T 6 


8.68 10" 5 


0.000274 


0.0455 


0.00204 


1.5 1(T 8 


1.29 10~ 6 


1.25 10~ 6 


0.00554 


3.11 lO" 5 



4.3 Simulations for e = 1 

When e = 0(1), the resolved and unresolved AP schemes are almost similar 
and we want to show that they give similar results as the resolved conven- 
tional scheme. In this way, we show that the AP scheme is as good as the 
conventional scheme when e = 0(1). We have already shown in the previous 
section that the former is much better than the latter when £< 1. 
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Figure 3: Comparison of the resolved conventional scheme (crosses), the 
resolved AP scheme (circles) and the exact drift-fluid limit (vertical bars) at 
t = 0.1 for e = 10~ 6 ; density n (top left), x-component of the momentum nu x 
(top right), y-component of the momentum nu y (bottom left), z-component 
of the momentum nu z (bottom right). 

Table 3: Maximum of relative difference between the computed solution and 
the exact drift-fluid limit (%) from the resolved AP scheme on the density 
n, x-component of the momentum nu x , ^-component of the momentum nu y , 
and maximum of absolute difference between the computed solution and the 
exact drift-fluid limit on the ^-component of the momentum nu z . 
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nu x 


nu y 


nu z 


10~ 5 


0.00873 


0.0074 


0.126 


0.017 


10~ 6 


8.72 10" 5 


0.000287 


0.0394 


0.00207 


1.5 10~ 8 


1.3 lO" 6 


1.25 10" 6 


0.0048 


3.16 10" 5 



When e — 1, the exact solution is no longer the drift-fluid limit solution. 
So, we restict ourselves to a comparison between the resolved AP scheme 
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Figure 4: Relative difference between the computed solution and the exact 
drift-fluid limit of the resolved AP scheme at time t = 0.1 for e = 10~ 6 ; den- 
sity n (top left), x-component of the momentum nu x (top right), y-component 
of the momentum nu y (bottom left). Absolute difference between the com- 
puted solution and the exact drift-fluid limit on the z-component of the 
momentum nu z (bottom right). 

Table 4: Maximum of relative difference between the computed solution and 
the exact drift-fluid limit (% ) from the non-resolved conventional scheme 
on the density n, x-component of the momentum nu x , y-component of the 
momentum nu y and maximum of absolute difference between the computed 
solution and the exact drift-fluid limit on the ^-component of the momentum 
nu z . 
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and the resolved conventional scheme. This comparison is shown in Fig. [9j 
We notice that the calculated solutions with the two schemes are indistin- 
guishable. Therefore for e of order of 1, both the resolved AP scheme and 
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Figure 5: Comparison of the non-resolved conventional scheme (crosses), 
the non-resolved AP scheme (circles) and the exact drift-fluid limit (vertical 
bars) at t = 0.1 for s = 10~~ 6 ; density n (top left), x-component of the 
momentum nu x (top right), y-component of the momentum nu y (bottom 
left), z-component of the momentum nu z (bottom right). 

Table 5: Maximum of relative difference between the computed solution and 
the exact drift-fluid limit (%) from the non-resolved AP scheme on the density 
n, x-component of the momentum nu x , y-component of the momentum nu y , 
and maximum of absolute difference between the computed solution and the 
exact drift-fluid limit on the ^-component of the momentum nu z . 
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resolved conventional scheme are comparable. 
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Figure 6: Relative difference between the computed solution and the exact 
drift-fluid limit of the non-resolved conventional scheme at time t = 0.1 for 
e = 10 -6 ; density n (top left), x-component of the momentum nu x (top 
right), y-component of the momentum nu y (bottom left). Absolute differ- 
ence between the computed solution and the exact drift-fluid limit on the 
^-component of the momentum nu z (bottom right). 

Table 6: Logarithms of the gyro-period r, maximum of time-steps used in 
the resolved AP scheme (AP) and non-resolved AP scheme (NAP) 
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AP 


NAP 


lir 5 
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-5.09 


-2.6 


lir 6 


-6 


-5.6 


-2.6 


1.5 1(T 8 


-7.83 


-6.51 
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4.4 Simulations for unprepared conditions 

For the sake of completeness, we show some numerical results obtained with 
unprepared boundary conditions. For this purpose, we introduce a second 
parameter e' = 10 -2 and use initial and boundary conditions as given by 
table [1] with e replaced by e' . On the other hand, e is kept at the value 
e = 10- 6 in the model fl2H), Q and in the scheme (EHJ), Q - 

On Fig. [TU] and [TTJ we display the values of the density and the three 
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Figure 7: Relative difference between the computed solution and the exact 
drift-fluid limit of the non-resolved AP scheme at time t — 0.1 for e = 
1CT 6 ; density n (top left), x-component of the momentum nu x (top right), y- 
component of the momentum nu y (bottom left). Absolute difference between 
the computed solution and the exact drift-fluid limit on the ^-component of 
the momentum nu z (bottom right). 

Table 7: CPU time (in s) used in the resolved conventional scheme (CONV) 
and non-resolved AP scheme (NAP) for computing the Euler-Lorentz model 
at final time t^ n (in s). Ratio of the CPU time of the conventional to the 
non-resolved AP schemes. 
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1(T 6 
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1584.21 


1.39 
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1.5 1(T 8 


0.01 


1149.54 


0.17 


6762 



components of the momentum as functions of x (resp. y) along the line 
y = 0.5 (resp. x = 0.5) for the resolved and non-resolved AP schemes and 
for the drift-fluid limit. The relative differences (absolute difference in the 
case of nu z ) of the solution with the drift-fluid limit are given on Fig. [T2lfor 
the resolved AP scheme and on Fig. [13] for the non-resolved one. 

Both the resolved and non-resolved AP schemes exhibit a significant dis- 
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Figure 8: Time-step (log scale) as a function of time for the resolved and 
non-resolved AP schemes when e = 10~ 6 . 

crepancy with the drift-fluid limit. This discrepancy originates in the ap- 
pearance of boundary layers which pollute the accuracy of the solution inside 
the domain. However, the discrepancy is much larger for the resolved AP 
scheme than for the non-resolved one. In many instances, the non-resolved 
AP scheme provides a fairly correct solution and its oscillations inside the 
boundary layers are less pronounced. This can be attributed to the larger 
time steps which provide a stronger relaxation rate towards the drift-fluid 
limit as well as a bigger amount of numerical diffusion than the small time- 
steps used in the resolved-AP schemes. 

These results show that the use of the non-resolved AP scheme in the 
case of non well-prepared boundary data at least provides a stable if not 
accurate solution. Additionally, it is expected that a suitable boundary layer 
analysis will permit to derive corrected boundary conditions which will take 
into account the influence of the boundary layer. The search for adequate 
boundary layer correctors will be the subject of future work. 



5 Conclusion 

The Euler-Lorentz model in the drift-fluid scaling for ion flow has been in- 
vestigated. First the drift-fluid limit has been studied and it has been shown 
that the parallel fluid velocity to the magnetic field is a solution of an el- 
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Figure 9: Comparison of the conventional scheme (crosses), and resolved AP 
scheme (circles) at t = 1 for 5 = 1; density n (top left), x-component of the 
momentum nu x (top right), y-component of the momentum nu y (bottom 
left), z-component of the momentum nu z (bottom right). 



liptic equation. Then, the full Euler-Lorentz model in the drift-fluid scaling 
has been investigated. A reformulation of the model has been provided, 
which shows that the parallel velocity is the solution of a wave equation with 
wave velocities tending to infinity as the scaling parameter e goes to zero. 
This reformulation allows to derive an Asymptotic Preserving scheme for the 
Euler-Lorentz model in the drift-fluid scaling, i.e. a scheme which is consis- 
tent with the full Euler-Lorentz model when e = 0(1) and which is consistent 
with its drift-fluid limit when e — > 0. The scheme allows to compute the so- 
lution of the Euler-Lorentz model when e< 1 with a time step independent 
of e. This property has been demonstrated numerically on a test example. It 
has been shown that the scheme is as good as the conventional scheme when 
e = 0(1) and that it provides a stable if not accurate solution in the case of 
non well-prepared boundary data. 

Forthcoming work will be devoted to the application of the AP scheme 
in the case of time varying and inhomogeneous magnetic fields, the coupling 
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Figure 10: Comparison of the resolved AP scheme (crosses), and unresolved 
AP scheme (circles) at t = 0.1 for unprepared initial and boundary condi- 
tions e = 10~ 6 and e' = 10~ 2 ; density n (top left), x-component of the 
momentum nu x (top right), y-component of the momentum nu y (bottom 
left), ^-component of the momentum nu z (bottom right). The computed so- 
lutions are shown for the section at middle y = 0.5 of the calculation domain 
Q along the z-direction. 



of the ion flow to the electron flow and to the electric and magnetic fields as 
well as to more rigorous stability analyses of the proposed schemes. 
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